using DataFrames
using QuadGK
using Distributions
using JLD
using Optim
using ForwardDiff
using CSV
using LinearAlgebra
using SpecialFunctions
using Random
using DelimitedFiles
using SparseArrays

clearconsole()

#-------------------------------------------------------------
# include Functions
#-------------------------------------------------------------
readDir = "$(pwd())/Functions/"
include(readDir *"vech.jl");
include(readDir *"VAR_Procedures.jl");
include(readDir *"Loaddata.jl");
include(readDir *"VAR_MoreLags_Procedures.jl")

#-------------------------------------------------------------
# choose specification files
#-------------------------------------------------------------
nDataSel  = "3"
nfVARSpec = "3"
nMDDSpec  = "1"
specDir   = "$(pwd())/SpecFiles/"
include(specDir * "/fVARspec" * nfVARSpec * ".jl")
include(specDir * "/MDDspec" * nMDDSpec * ".jl")


#-------------------------------------------------------------
# load aggregate data
#-------------------------------------------------------------
dataDir  = "$(pwd())/data/Para" * nDataSel *"/"
agg_data = loadaggdata(1,Tend)
n_agg    = size(agg_data)[2]

#-------------------------------------------------------------
# Iterate over number of knots in specification file
#-------------------------------------------------------------
sNameLoadDir = "fVAR" * nfVARSpec
sNameSaveDir = "fVAR" * nfVARSpec * "_MDD" * nMDDSpec
loaddir  = "$(pwd())/results/Para" *nDataSel * "/" * sNameLoadDir *"/";
savedir  = "$(pwd())/results/Para" *nDataSel * "/" * sNameSaveDir *"/";
try mkdir(savedir) catch; end

MDD_Laplace_sum         = zeros(hyper_n, K_vec_n+size(hyper_grid)[2])
Pen_Laplace_sum         = zeros(hyper_n, K_vec_n+size(hyper_grid)[2])
lambda_hat_K_LaplaceMax = zeros(K_vec_n,7)

for ii = 1:K_vec_n

    time_init_loop = time_ns()
    K = K_vec[ii]
    print("Approximation order K $K \n")

    #-------------------------------------------------------------
    # Load coefficients from density estimation
    #-------------------------------------------------------------

    PhatDensCoef_factor, MDD_GoF, VinvLam_all, ~ , ~  = loaddensdata(1,Tend,K,nfVARSpec)
    n_cross = size(PhatDensCoef_factor)[2]

    #-------------------------------------------------------------
    # Compute MDD term2 over hyperparameter grid
    #-------------------------------------------------------------
    Ktilde      = size(PhatDensCoef_factor)[2]

    Pen_Laplace = zeros(Tend-Tdrop,1)
    for tt = 1:(Tend-Tdrop)
        Pen_Laplace[tt,1] = -0.5*logdet(VinvLam_all[:,:,Tdrop+tt]) + Ktilde/2*log(2*pi)
    end

    for ll=1:hyper_n

        lambda1 = hyper_grid[ll,1]
        lambda2 = hyper_grid[ll,2]
        lambda3 = hyper_grid[ll,3]
        nlags   = Int(hyper_grid[ll,4])

        # adjust the sample
        agg_data_adj            = agg_data[Tdrop-nlags+1:end,:]
        PhatDensCoef_factor_adj = PhatDensCoef_factor[Tdrop-nlags+1:end,:]
        MDD_GoF_est             = MDD_GoF[Tdrop+1:Tend,:]

        # Combine aggregate and cross sectional observations and create YY XX for VAR estimation
        data_adj = [ agg_data_adj PhatDensCoef_factor_adj ]
        YY, XX, PHIols, Sols, SIGMAols = data_to_YY(data_adj, nlags, nlags)

        # generate prior
        igammadf = size(YY)[2]+sigdf
        lambda  = [lambda1 lambda2 lambda3 2]
        prior    = prior_ACP_stru(YY, SIGMAols, n_agg, igammadf, nlags, lambda)

        # Define some constants
        n       = n_agg+n_cross
        p       = nlags

        if ii == 1
            MDD_Laplace_sum[ll,1:4] = hyper_grid[ll,:]'
            Pen_Laplace_sum[ll,1:4] = hyper_grid[ll,:]'
        end

        lnpYY = ml_VAR_ACP(nlags,YY,XX,prior)

        Pen_Laplace_sum[ll,4+ii] = sum(Pen_Laplace[:,1]) + lnpYY
        MDD_Laplace_sum[ll,4+ii] = sum(MDD_GoF_est) + Pen_Laplace_sum[ll,4+ii]

    end

    lambda_hat_K_LaplaceMax[ii,:] = [ K n_cross MDD_Laplace_sum[argmax(MDD_Laplace_sum[:,4+ii]),[1 2 3 4 (4+ii)]] ]

    time_loop=signed(time_ns()-time_init_loop)/1000000000
    print("Elapsed time           = $(time_loop) \n")
    print("================================= \n")

end

sNameFile = "fVAR" * nfVARSpec * "_MDD" * nMDDSpec
CSV.write(savedir * sNameFile * "_MDD_Laplace_sum.csv", DataFrame(MDD_Laplace_sum,:auto))
CSV.write(savedir * sNameFile * "_Pen_Laplace_sum.csv", DataFrame(Pen_Laplace_sum,:auto))

sNameCols = ["K", "Ktilde", "Lam1", "Lam2", "Lam3", "Nlags", "MDD"]
CSV.write(savedir * sNameFile * "_lambda_hat_K_LaplaceMax.csv", DataFrame(permutedims(sNameCols), :auto), header=false)
CSV.write(savedir * sNameFile * "_lambda_hat_K_LaplaceMax.csv", DataFrame(lambda_hat_K_LaplaceMax, :auto),append=true)
